{
 "nbformat": 4,
 "nbformat_minor": 2,
 "metadata": {
  "language_info": {
   "name": "python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "version": "3.7.6-final"
  },
  "orig_nbformat": 2,
  "file_extension": ".py",
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3,
  "kernelspec": {
   "name": "python37664bit8b61c5f9b63a430188c2e600c83aebf6",
   "display_name": "Python 3.7.6 64-bit"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "二分法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisection(f, a, b, eps, eta=1e-16, verbose=False):\n",
    "    \"\"\"二分法求根\n",
    "    \n",
    "    对方程 f(x) = 0 在区间 [a, b]，使用二分法求根。\n",
    "    做 ceil((log(((b - a) / eps), 2) - 1) 次迭代，使结果满足精度 eps。\n",
    "    \n",
    "    实现参考: 数值分析[谷根代，杨晓忠 等著]2011年版.P18.算法1\n",
    "\n",
    "    Args:\n",
    "        f: function, 一元函数，表示要求根的方程: f(x) = 0.\n",
    "        a, b: float, 有根区间 [a, b] 的端点.\n",
    "        eps: float, 给定精度.\n",
    "        \n",
    "        eta: float, 当 abs(f(x)) <= eta 时停止计算, default 1e-16.\n",
    "        verbose: bool, 打印出二分法计算的表格, default False.\n",
    "\n",
    "    Returns:\n",
    "        (x_final, N)\n",
    "\n",
    "        x_final: float, 二分法求得的近似根\n",
    "        N: int, 迭代次数\n",
    "\n",
    "    Raises:\n",
    "        ValueError: 给定的 f(a) * f(b) < 0 时无法求解，抛出异常\n",
    "    \"\"\"\n",
    "    if f(a) * f(b) > 0:\n",
    "        raise ValueError(\"rootless interval: f(a) * f(b) < 0\")\n",
    "\n",
    "    N = math.ceil((log(((b - a) / eps), 2) - 1).evalf(5))\n",
    "\n",
    "    if verbose:\n",
    "        print(f'N = {N}\\n')\n",
    "    \n",
    "    x = (a + b) / 2\n",
    "    \n",
    "    if verbose:\n",
    "        print(f'n \\t (a, b) \\t f(x_n)')\n",
    "        print('-'*35)\n",
    "    \n",
    "    for n in range(N+1):\n",
    "        if abs(f(x)) <= eta:  # f(x) == 0\n",
    "            break\n",
    "\n",
    "        if verbose:\n",
    "            print(f'{n} \\t ({a}, {b}) \\t f({x})={f(x)}')\n",
    "        \n",
    "        if f(x) * f(a) < 0:\n",
    "            b = x\n",
    "        else:\n",
    "            a = x\n",
    "        x = (a + b) / 2\n",
    "        n += 1\n",
    "        \n",
    "    if verbose:\n",
    "        print(f'{n} \\t ({a}, {b}) \\t -')\n",
    "    \n",
    "    x_final = (a + b) / 2\n",
    "    if verbose:\n",
    "        print(f'\\nresult: x = ({a}+{b})/2 = {x_final}')\n",
    "\n",
    "    return x_final, N+1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "N = 10\n\nn \t (a, b) \t f(x_n)\n-----------------------------------\n0 \t (0, 1) \t f(0.5)=0.733635780821064\n1 \t (0.5, 1) \t f(0.75)=0.263094345458695\n2 \t (0.75, 1) \t f(0.875)=0.0661805371209897\n3 \t (0.875, 1) \t f(0.9375)=-0.0228698549070950\n4 \t (0.875, 0.9375) \t f(0.90625)=0.0208763999947352\n5 \t (0.90625, 0.9375) \t f(0.921875)=-0.00119109771321235\n6 \t (0.90625, 0.921875) \t f(0.9140625)=0.00979401298210625\n7 \t (0.9140625, 0.921875) \t f(0.91796875)=0.00428930379438952\n8 \t (0.91796875, 0.921875) \t f(0.919921875)=0.00154606529293533\n9 \t (0.919921875, 0.921875) \t f(0.9208984375)=0.000176724441991016\n10 \t (0.9208984375, 0.921875) \t f(0.92138671875)=-0.000507376461447939\n11 \t (0.9208984375, 0.92138671875) \t -\n\nresult: x = (0.9208984375+0.92138671875)/2 = 0.921142578125\n"
    },
    {
     "data": {
      "text/plain": "(0.921142578125, 11)"
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return 2 * exp(-x) - sin(x)\n",
    "\n",
    "bisection(f, 0, 1, 0.0005, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dichotomy(f, a, b, eps, eta=1e-16, verbose=False):\n",
    "    \"\"\"二分法求根\n",
    "\n",
    "    对方程 f(x) = 0 在区间 [a, b]，使用二分法求根。\n",
    "    一直迭代到 abs(b - a) <= eps 为止。\n",
    "\n",
    "    实现参考：《实验二  非线性方程求根》\n",
    "\n",
    "    Args:\n",
    "        f: function, 一元函数，表示要求根的方程: f(x) = 0.\n",
    "        a, b: float, 有根区间 [a, b] 的端点.\n",
    "        eps: float, 根的容许误差.\n",
    "        \n",
    "        eta: float, abs(f(x)) 的容许误差, default 1e-16.\n",
    "        verbose: bool, 打印出二分法计算的表格, default False.\n",
    "\n",
    "    Returns:\n",
    "        (x_final, N)\n",
    "\n",
    "        x_final: float, 二分法求得的近似根\n",
    "        N: int, 迭代次数\n",
    "\n",
    "    Raises:\n",
    "        ValueError: 给定的 f(a) * f(b) < 0 时无法求解，抛出异常\n",
    "    \"\"\"\n",
    "    if f(a) * f(b) > 0:\n",
    "        raise ValueError(\"rootless interval: f(a) * f(b) < 0\")\n",
    "    \n",
    "    if verbose:\n",
    "        print(f'n \\t (a, b) \\t f(x_n)')\n",
    "        print('-'*35)\n",
    "    \n",
    "    n = 0\n",
    "    while abs(b - a) > eps:\n",
    "        x = (a + b) / 2\n",
    "        fx = f(x)\n",
    "        if verbose:\n",
    "            print(f\"{n}\\t {a, b}\\t f({x})={fx}\")\n",
    "\n",
    "        if abs(fx) <= eta:\n",
    "            break\n",
    "        \n",
    "        if f(a) * fx < 0:\n",
    "            b = x\n",
    "        else:\n",
    "            a = x\n",
    "        \n",
    "        n += 1\n",
    "\n",
    "    x_final = (a + b) / 2\n",
    "    if verbose:\n",
    "        print(f\"{n}\\t {a, b}\\t -\")\n",
    "        print(f'\\nresult: x = ({a}+{b})/2 = {x_final}')\n",
    "    \n",
    "    return x_final, n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "(0.921142578125, 11)"
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dichotomy(f, 0, 1, 0.0005, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fixed_point_iter(phi, x_0, max_steps=25, verbose=False):\n",
    "    \"\"\"不动点迭代\n",
    "\n",
    "    Args: \n",
    "        phi: function, 迭代函数\n",
    "        x_0: float, 初值\n",
    "        \n",
    "        max_steps: int, 最大迭代次数\n",
    "        verbose: bool, 打印出每一步的值，default False.\n",
    "\n",
    "    Returns: \n",
    "        x_final: float, 最终的近似根 x \n",
    "    \"\"\"\n",
    "    x = x_0\n",
    "    for i in range(max_steps):\n",
    "        if verbose:\n",
    "            print(f'x_{i} \\t {x}')\n",
    "        x = phi(x)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "x_0 \t 1.5\nx_1 \t 1.4444444444444444\nx_2 \t 1.4792899408284024\nx_3 \t 1.456976\nx_4 \t 1.4710805833200253\nx_5 \t 1.4620905354712408\nx_6 \t 1.4677905760195855\nx_7 \t 1.464164380462178\nx_8 \t 1.4664663557170745\nx_9 \t 1.465003040566855\n"
    },
    {
     "data": {
      "text/plain": "1.4659324390818347"
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fixed_point_iter(lambda x: 1 + 1 / x**2, 1.5, max_steps=10, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton_iter(f, x_0, eps=0, eta=0, df=None, max_steps=20, frac=False, verbose=False):\n",
    "    \"\"\"牛顿迭代法\n",
    "\n",
    "    Args: \n",
    "        f: function, 迭代函数\n",
    "        x_0: float, 初值\n",
    "        eps: float, 根的容许误差, default 0.\n",
    "        eta: float, abs(f(x)) 的容许误差, default 0.\n",
    "        \n",
    "        df: function, f 的导函数, 默认 None 表示自动调用 sympy.diff 求导(会导致后续迭代中使用分数运算)。\n",
    "        max_steps: int, 最大迭代次数, default 20.\n",
    "        frac: bool, True 则输出分数(仅对 df=None 时生效)，否则使用 float, default False.\n",
    "        verbose: bool, 打印出每一步的值, default False.\n",
    "\n",
    "    Returns: \n",
    "        x_final: float, 最终的近似根 x\n",
    "    \"\"\"\n",
    "\n",
    "    if df == None:\n",
    "        __x = Symbol('x')\n",
    "        __df = diff(f(__x), __x)\n",
    "        df = lambda x: __df.subs(__x, x)\n",
    "\n",
    "    x = x_0\n",
    "    for i in range(max_steps):\n",
    "        if verbose:\n",
    "            print(i, x)\n",
    "        \n",
    "        x_next = x - f(x) / df(x)\n",
    "        \n",
    "        if abs(x_next - x) <= eps or abs(f(x)) <= eta:\n",
    "           break\n",
    "        \n",
    "        x = x_next\n",
    "    \n",
    "    if not frac:\n",
    "        x = float(x)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "0 2\n1 2.1\n2 2.094568121104185\n"
    },
    {
     "data": {
      "text/plain": "2.094568121104185"
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return x ** 3 - 2 * x - 5\n",
    "\n",
    "def df(x):\n",
    "    return 3 * x ** 2 - 2\n",
    "\n",
    "# newton_iter(f, 2, eps=0.5e-3, frac=True, verbose=True)\n",
    "newton_iter(f, 2, eps=0.5e-3, df=df, frac=True, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "def single_point_truncation(f, x_0, x_1, eps=0, eta=0, max_steps=20, verbose=False):\n",
    "    \"\"\"单点弦截法\n",
    "\n",
    "    Args: \n",
    "        f: function, 迭代函数\n",
    "        x_0, x_1: float, 初值\n",
    "        eps: float, 根的容许误差, default 0.\n",
    "        eta: float, abs(f(x)) 的容许误差, default 0.\n",
    "        \n",
    "        max_steps: int, 最大迭代次数, default 20.\n",
    "        verbose: bool, 打印出每一步的值, default False.\n",
    "\n",
    "    Returns: \n",
    "        x_final: float, 最终的近似根 x\n",
    "    \"\"\"\n",
    "\n",
    "    f_x0 = f(x_0)\n",
    "    \n",
    "    x = x_1\n",
    "    for i in range(max_steps):\n",
    "        if verbose:\n",
    "            print(i, x)\n",
    "        \n",
    "        x_next = x - f(x) / (f(x) - f_x0) * (x - x_0)\n",
    "        \n",
    "        if abs(x_next - x) <= eps or abs(f(x)) <= eta:\n",
    "           break\n",
    "        \n",
    "        x = x_next\n",
    "    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "0 1\n1 2.2\n2 2.088967971530249\n3 2.094861151990966\n"
    },
    {
     "data": {
      "text/plain": "2.094861151990966"
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return x ** 3 - 2 * x - 5\n",
    "\n",
    "single_point_truncation(f, 2, 1, eps=0.5e-3, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "def secant_method(f, x0, x1, eps=0, eta=0, max_steps=20, verbose=False):\n",
    "    \"\"\"两点弦截法（割线法）\n",
    "\n",
    "    Args: \n",
    "        f: function, 迭代函数\n",
    "        x0, x1: float, 初值\n",
    "        eps: float, 根的容许误差, default 0.\n",
    "        eta: float, abs(f(x)) 的容许误差, default 0.\n",
    "        \n",
    "        max_steps: int, 最大迭代次数, default 20.\n",
    "        verbose: bool, 打印出每一步的值, default False.\n",
    "\n",
    "    Returns: \n",
    "        x_final: float, 最终的近似根 x\n",
    "    \"\"\"\n",
    "    for i in range(max_steps):\n",
    "        if verbose:\n",
    "            print(i, x1)\n",
    "        \n",
    "        x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))\n",
    "        x0, x1 = x1, x2\n",
    "\n",
    "        if abs(x1 - x0) <= eps or abs(f(x1)) <= eta:\n",
    "           break\n",
    "            \n",
    "    return x1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "0 30\n1 22.8\n2 24.545454545454547\n3 24.746543778801843\n4 24.73860275369709\n"
    },
    {
     "data": {
      "text/plain": "24.738633748750722"
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "secant_method(lambda x: x ** 2 - 612, 10, 30, max_steps=5, verbose=True) # Root: 24.738633748750722"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "题目 \n",
    "求方程 $f(x)=x^3+x^2-3x-3=0$ 在 1.5 附近的根.（误差限为 $\\epsilon=1e-6$, $\\eta=1e-9$）\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": "牛顿迭代法:\n0 1.5\n1 1.77777777777778\n2 1.73336066694000\n3 1.73205192940947\n4 1.73205080756970\n单点弦截法:\n0 2\n1 1.6923076923076923\n2 1.7390156515180086\n3 1.7308625826467308\n4 1.7322544663053168\n5 1.7320159286895187\n6 1.7320567817872679\n7 1.7320497843005194\n8 1.732050982835706\n两点弦截法:\n0 2\n1 1.6923076923076923\n2 1.7257977285018928\n3 1.7322172842612025\n4 1.732050123979108\n\nresult:\n牛顿迭代法: 1.7320508075697012\n单点弦截法: 1.732050982835706\n两点弦截法: 1.7320508074943775\n"
    }
   ],
   "source": [
    "def f(x):\n",
    "    return x ** 3 + x ** 2 - 3 * x - 3\n",
    "\n",
    "x0 = 1.5\n",
    "x1 = 2\n",
    "eps = 1e-6\n",
    "eta = 1e-9\n",
    "\n",
    "result = {}\n",
    "print(\"牛顿迭代法:\")\n",
    "result[\"牛顿迭代法\"] = newton_iter(f, x0, eps=eps, eta=eta, verbose=True)\n",
    "print(\"单点弦截法:\")\n",
    "result[\"单点弦截法\"] = single_point_truncation(f, x0, x1, eps=eps, eta=eta, verbose=True)\n",
    "print(\"两点弦截法:\")\n",
    "result[\"两点弦截法\"] = secant_method(f, x0, x1, eps=eps, eta=eta, verbose=True)\n",
    "\n",
    "print(\"\\nresult:\")\n",
    "for k in result:\n",
    "    print(f'{k}: {result[k]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}